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ABSTRACT 

Numerical  methods  are  described  for  the  calculation 
of  two-dimensional  time  dependent  inviscid  flows  which 
contain  shocks.   The  problems  were  considered  to  be  time 
dependent  and  their  solutions  are  found  asymptotically 
for  large  values  of  time   (t  ->  oo ) .   Results  are  given 
for  oblique  shock  reflection  and  Mach  reflection.   The 
numerical  scheme  was  tested  for  values  of  Mach  numbers 

that  ranged  between  1.3  -^  M   <  3.0. 

—  00  —  ^ 

The  conservation  equations  describing  such  flows 
are  integrated  over  a  region  which  contains  multiple  shocks 
Two  of  the  difference  methods  described  are  high  order 
accuracy  schem.es  (the  method  of  Lax  and  Wendroff)  which 
closely  approximate  the  jump  conditions  across  disconti- 
nuities and  gives  the  shock  transition  within  three  mesh 
widths.   The  initial  value  problems  that  were  considered 
had  initial  data  close  to  the  asymptotic  solution  and  data 
that  were  arbitrarily  rem.ote  from  the  desired  solution. 
The  methods  used  were  able  to  handle  the  severe  transients 
encountered  in  all  cases. 
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In  addition  to  the  above  methods  two  other  difference 
analogs  to  the  partial  differential  equations  of  motion 
were  tested  and  compared.   Numerical  results  are  given  for 
the  oblique  reflection  and  Mach  reflection  calculation. 
Comparison  of  the  computed  data  with  the  exact  solution  of 
ordinary  reflection  is  good.   The  behavior  of  the  Mach 
reflection  calculation  seems  to  be  qualitatively  correct. 


-      3 


TABLE  OP  CONTENTS 

Abstract  2 
Section 

1.  Introduction  $ 

2.  The  Differential  Equations  Y 
3-   Lax-Wendroff  Difference  Equations  ^ 

4.  A  Third  Difference  Scheme  of 

Second  Order  Accuracy  15 

5.  A  Fourth  Difference  Scheme  Using 

Pseudo  Viscosity  l6 

6.  Regular  Net  Points  19 

7.  Boundary  Points  20 

8.  Shock  Relations  24 

9.  Initial  Conditions  28 

10.  Numerical  Tests  of  Stability  29 

11.  Numerical  Results  3I 
Bibliography  33 


-  4 


NYO-10,433 

NUMERICAL  CALCULATIONS  OF  MULTIDIMENSIONAL 

SHOCKED  FLOWS* 

by 

Samuel  Z.  Bursteln 

1 .   Introduction 

With  the  advent  of  high  speed  computing  machines,  It 
has  become  feasible  to  generate  detailed  solutions  to 
multidimensional  fluid  dynamics  problems.   Some  of  these 
problems  deal  with  flows  containing  discontinuities  called 
shocks.   Shocks  may  be  formed  in  the  fluid  in  many  ways, 
e.g.,  spontaneously  as  a  natural  result  of  the  nonlinear 
fluid  motion,  as  a  result  of  explosions  or  combustion 
processes,  or  as  a  result  of  bodies  causing  an  obstruction 
or  restriction  in  the  flow.   This  paper  will  be  concerned 
with  those  flows  containing  shocks  formed  as  a  result  of 
the  last  category.   In  particular,  the  cases  or  ordinary 
oblique  shock  reflection  and  that  of  Mach  reflections  are 
considered.   The  flows  considered  here  are  steady  in  time. 
They  are  generated  by  assuming  that  the  solutions  are  time 
dependent  and  determining  asymptotic  solutions  as   t  ->  oo  . 
Hence,  the  solutions  contain  the  dynamics  of  the  approach 


*The  work  presented  in  this  paper  is  supported  by  the 
AEC  Computing  and  Applied  Mathematics  Center,  Courant 
Institute  of  Mathematical  Sciences,  New  York  University, 
under  Contract  AT( 30-1) -l480  with  the  U.S.  Atomic  Energy 
Commission. 
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to  the  asymptotic  solution  as  well  as  the  steady  state. 

Until  recently,  in  most  computations,  the  equations 
of  motion  were  differenced  from  the  Eulerlan  or  Lagranglan 
form.   Peter  Lax  [1-4]  suggested  that  the  conservation  or 
divergence  form  of  the  partial  differential  equations 
be  used  to  generate  the  difference  equations.   It  Is 
with  these  equations,  and  some  modified  difference 
equations,  that  numerical  results  will  be  given  for  the 
two  dimensional  cases  considered.   In  all,  four  finite 
difference  analogs  were  tried.   Three  are  second  order 
schemes  while  the  fourth  Is  a  first  order  method  using  a 
pseudo  viscous  pressure  for  damping. 

This  work  has  been  motivated  by  Professor 
Robert  D.  Rlchtmyer,  who  has  presented  a  program  of  fluid 
dynamics  research,  which  has  as  one  of  its  goals  the  more 
successful  use  of  high  speed  computing  machines  In  fluid 
dynamic  research  [5].   The  author  would  like  to  thank 
John  Gary  and  Professor  Herbert  Keller  for  their  stimulate 
Ing  discussions,  and  Professor  Peter  Lax  for  many  helpful 
ideas,  including  the  suggestion  of  Gilbert  Strang. 
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2.   The  Differential  Equations 

The  Invlscld  equation  of  fluid  motion  may  be  written 
In  divergence  form,  fulfilling  the  principle  of  conserva- 
tion of  mass,  momentum  and  energy: 


pt  =  -^x  -  ^y  ' 

2 
f      ,  m  >     /nm* 

-t    =    -^P  +T^--  ^T^y    ' 


(2.1) 


,t     '  p  ,x   '^  p  ,y 

h      =      -[^(P  +  E)]    -  [^(p+E)] 


y 


Here  p,   m  =  pu,   n  =  pv  are  the  mass  per  unit  volume 

and  momenta  per  unit  volume  In  the  x  and  y  directions, 

12    2 
E  =  p[e  +-^(u  +  V  )]   Is  the  total  energy  per  unit  volume. 

The  subscripts  following  the  comma  denote  partial  differ- 
entiation.  The  pressure  p,   the  fifth  unknown,  may  be 
expressed  in  terms  of  the  density  and  specific  internal 
energy,  by  the  equation  of  state,   p  =  P(p,e).   In  this 
paper,  the  functional  P   is  given  by  the  ideal  gas  law 

p   .   p(^_i)e  =   (^_i)[E-linLJLn!)] 

so  that  the  conservation  equations  (2.1)  may  be  written 
completely  in  terms  of  the  conservation  variables  p,  m, 
n,  E. 

If  w  is  a  vector  quantity  which  is  defined  in  some 
volume  region  of  space,  its  time  rate  of  change  depends 
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upon  the  flux,  represented  by  the  vectors   f  and  g, 
into  this  space.   This  statement  leads  to  the  conserva- 
tion law  in  divergence  form,  i.e., 


(2.2) 


^t  -   ^x+S,y  ' 


where  the  vector  function  w  has  the  four  components 


w  = 


and  the  given  vectors   f  and  g  have  the  form 


f(w)   = 


-  m 


7-3  1^ 
2   P 


Em 
^  P 


(r-DE  +  ^ljii^ 


mn 

P 

3     2 
7-1  m-^  +  mn 


g(w)   = 


-  n 

nm 
P 

2   p     ^^  ^^^  ^   2   p 

3     2 
En   7-1  n-^+  nm 

~  ^  p  "   2      2 


It  is  this  differential  equation  (2.2)  and  its  ability  to 
predict  so-called  weak  solutions,  solutions  which  contain 
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shocks,  which  Is  the  object  of  this  paper.   A  classical 
solution  to  a  fluid  dynamic  problem  is  a  solution  which 
consists  of  the  description  of  w,   the  vector  valued 
function,  up  to  a  time  when  shocks  appear  in  the  domain 
of  interest-   A  weak  solution  is  a  solution  which  may  be 
continued  so  as  to  predict  the  shocked  region  for  all 
subsequent  times.   The  continuation  of  the  solution  is 
possible  since  the  Rankine-Hugoniot  condition  is  satisfied 
across  the  discontinuities  predicted  by  the  weak  solution 
(see  [3]) . 

The  point  is  that  the  approach  taken  here  does  not 
consider  the  moving  shock  as  an  interior  boundary  (i.e., 
the  approach  is  related  to  that  taken  in  [6])  as  compared 
with  the  method  of  treatment  of  R .  D.  Richtmyer  [?]. 

3.   Lax-Wendroff  Difference  Equations 

The  Lax-Wendroff  method  [3-4]  is  based  on  the  Taylor 

expansion  of  the  vector  function  w(x,y,t+At)   so  as  to 

1   2 
include  the  second  order  term:   -?;  At  w  ^,  .   This 

d  ,  ut 

accomplishes  two  objectives: 

A)  It  is  this  term  which,  when  added  to  the  "intuitive" 

differencing  scheme  (based  on  equation  (2.2)) 

(3.1)     wCx,y,t+At)   =  w(x,y,t)  +  [f,.(t)+  g,(t)]At 

X     y 

produces  a  stable  difference  analog.   It  is  a  well-known 
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fact  that  the  above  differencing  scheme  (3>1)  is 
unconditionally  unstable.   The  central  difference  approxi- 
mations could  however  be  replaced  by  a  forward  or  backward 
difference,  resulting  in  a  difference  scheme  which  is 
conditionally  stable,  but  only  accurate  to  0(At  ). 

f/s  =   „    .    ( f  ( x+Z\x )    -   f  (x-Ax) )      denotes   the   centered 
X        2  Ax      ^  ^ 

difference  quotient  which  approximates  the  derivative  with 

a  truncation  error  which  is   0(Ax  ) . 

B)  It  increases  the  order  of  the  truncation  error  to 

"Z 

OCAt"^)  .   Lax  and  Wendroff  observed  that  the  form  of  this 
added  second  order  term  may  be  established  by  the 
differential  equation  itself 

^^t\tt   =  I  ^t'(^x-^^y^t   =  l^t'ff,tx+^ty) 
(3.2)  =  I  At2[(Aw  ^)  ^  +  (Bw  , )  ^1 

=  \   At2[(A(f   +  g  ^))  ^  +  (B(f   +  g  „))  ^] 

Here  we  have  used  A  and  B   to  represent  the  Jacoblans 
of   f(w)  and  g(w)   with  respect  to  w. 

Hence,  we  can  write  the  unknown  vector  in  terms  of 
the  previously  known  spatial  variations 


w(x,y,t+At)   =   w(x,y,t)  +  (f^+  g^^)At 

(5.5)  +  |[A(f^+gy)]^+  [B(f^+gy)]y| 

+  O(At^)  . 
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At2 


Lax  and  Wendroff  modified  this  relation  by  adding  a  term 
0(At  )   which  for  smooth  flows  Increases  the  allowable  time 
Increment.   This  additional  stabilizing  term  is 

i3'^)  -   ^  At^CA^  +  B^)w 

o    ^        ,xxyj 

Equation  (3 •3)  and  equations  (3«3)  plus  (3.4)  give  two 
difference  methods  for  the  solution  of  equation  (2.2). 
They  are  nine-point  difference  schemes  since  they  use  the 
point   (x,y)  and  its  eight  nearest  neighbors.   Although 
there  is  no  unique  way  of  differencing  equation  (3.3),  in 
Section  6  we  show  the  conservation  form  of  differencing. 

The  linear  stability  limit  for  (3.3)  may  be  given  in 
terms  of  the  time  and  space  Increment  ratios  and  the  flow 

u  +v  /c   (c   is  the  local  sound  speed) 

(3.5)      c^  <  smrl ,      ef  <  smrl 

Ax     -         ^  Ay     -         ^ 

and  the  condition  for  (3-3)  plus  (3-4)  when  Zix  =  Ay  =  A 
is 

(3-6)     ^     -    yS   '  ^  =  eigenvalue   JA^  +  B^]   . 

max      ^       / 
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4.   A  Third  Difference  Scheme  of  Second  Order  Accuracy 

The  amplification  matrix  G-,   (see  [9])  of  the  last 
two  difference  schemes  can  be  given  by 

(4.1)      G^(^,ii)   =   I+K^A  +  TiB)  -  |(^A+iiB)^  +  O(l^) 


and 


G(|,Ti)   =   G^(^,ii)  -lik^+B^)i%^ 


This  was  a  result  of  expanding  the  amplification  matrix 

of  exact  solutions  of  (2.2),  Q^^i^-^^i    ^      for  small  | 

and  ri   which  are  the  wave  numbers  equal  to  2ir/0    , 

2ir/6    ,      where  0        and  B        are  the  wave  lengths  In  the 
^  y  X        y 

X  and  y  directions.   Hence,  In  order  to  generate 
additional  difference  schemes  (possibly  with  Increased 
stability),  one  may  try  to  find  an  approximation  to  the 
exact  amplification  matrix.   The  approximation  should 
differ  from  the  exact  G   only  in  third  or  higher  order 
terms.   This  Idea  was  suggested  by  G.  Strang.   It  Is  not 
difficult  to  show  that 

(4.3)     giUA+nJjj   ^  e e +^ e ^  q^^^)^   ^ 

Before  expanding  the  right-hand  side  of  (4.^),  allow  each 
of  the  four  elemental  matrices  to  be  expanded  In  a  Taylor 
series.   This  procedure  yields  for  the  right-hand  side 
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(kA)  I  +l(4A4^B)  -  -I  i^A^   -  ^^   iri    -  \   Ti^B^  +  k 


with 

„      A^'b  +B  A   ^2^2    .  ,AB''+  B  A  t^2  ^  A^+BA   .2, 

The  following  approximations  are  used  correct  to  0(^'^): 


(4.5) 


2 
^   =   sin  ^  ,     ^    =   2(l-cos  ^) ,  ^ri  =  sln|  slnr), 

2 

■q   =   sin   T)  ^         1    =   2(1  -  cos  ■q )  . 


They  are  substituted  into  (4.4)  and  the  resulting  function 
is  denoted  by  Gp(C,Ti)j  the  amplification  matrix  of  a  new 
difference  scheme: 

Gg  =  I  +  i(A  sin  e  +  B  sin  r\)    -   A^(l  -  cos  i) 

-  B^(l   -    cos  T])    -  ^g^^   sin  i    sin  t^ 
(^•^^         ,  aS^+bV  ,,  ,,    ,,  , 

+  ^ (1-cos^)  (1-cosri) 

.  ,A^B+BA^   .    ,-,      .,  ,  AB^+B^A  .       .r■^  w 

-  i( 2 s^^  ri(l-cos  ^)  H ^ sm  ^(1-cos  t]  )  ) 

Equation  (4.3)  and  hence  equation  (4.6)  differ  from  the 
amplification  matrix  associated  with  the  first  two 
difference  schemes  only  by  terms   0(^^).   Clearly,  the 
difference  scheme  associated  with  equation  (4.6)  is  also 
accurate  to  Q(4^).   The  additional  terms  are  the 
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difference  analogs  of  the  differential  operators 

^(AB^+B^A)  ^  +^(a2b  +  Ba2)  1— 

hxby  hx   Sy 

Sx  Sy 

and  will  be  referred  to  as  the  triple  viscosity,  a 
stabilizing  term.   This  result  can  be  seen  If  (4.6') 
operates  on  a  solution  of  the  form   e  ^^   ^"^  .   The 
value  of  the  three  additional  terms  lies  In  the  fact  that 
the  amplification  matrix  may  be  written  In  the  separated 
form 

(4.7)  G,(e,,)   ==   M(ON(T^)  +N(T^)M(0 


where  we  use  the  abbreviations 


(4.8)  M(0   =   e^^^  ,  N(ti)   =   e^"^^ 


It  is  a  well-known  fact  of  hyperbolic  equations  that  A 
and  B  are  real,  and  Prledrlchs  has  shown  that  A  and 
B   can  be  symmetrized.   Then  M=M^=M,   I.e.,   M  is 
equal  to  its  transposed  complex  conjugate  M   as  a 
result  of  the  following  relations 

2.2 


M  =  I  +  1  A|  -  ^-^  +  Oii^) 
M*  =   I  -  lA^I  -  (A^)'^  -^  +  0(1^: 
=   I  -  1A|  -  AJ_  +  0(4^) 
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which  means  that  M  is  Hermltlan  and  also  MM  =  I  =  M  M, 
I.e.,  the  matrices  denoted  by  M   (and  N)   are  unitary. 
A  unitary  transformation  corresponds  to  a  generalized 
rotation  which  leaves  all  vectors  unaltered  in  length. 
We  can  now  form  the  inner  product  using  equation 
(4.7).   If  q  is  an  arbitrary  unit  vector,  then 

2(G2(e,Ti)q,q)   =   ( (M(e  )N(ti  )  +  N(ii  )M(e  )  )q,q)  . 

The  unitary  property  says 

(M(e)q,q)  <      1      . 

A  similar  statement  can  be  made  for  N('n).   Hence,  the 
amplification  matrix  (4.7)  satisfies  the  inequality 

According  to  Lax's  new  stability  theorem  [4],  this  implies 
that  the  absolute  value  of  all  integer  powers  of  the 
amplification  matrix  is  bounded,  and  as  a  result,  the 
associated  difference  equations  are  stable  if  the  time 
increment  satisfies  the  eigenvalue  Inequalities 

(^.10)  |A(A)|  ^      <      1    ,  |A(B)|  ^   <   1   , 

max  max   "^ 

which  is  the  Courant-Friedrichs-Lewy  condition  (8).   It 
is  clear  then  that  the  three  additional  terms  in  (4.6) 
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act  In  a  stabilizing  manner  so  that  the  allowable  step 
(4.10)  may  be  maximiim.   The  price  one  pays,  though,  is 
the  increase  in  computation  time  per  net  point  to  the 
matrix  multiplications. 

5.   A  Fourth  Difference  Scheme  Using  Pseudo  Viscosity 

The  idea  of  introducing  an  artificial  viscosity 
into  shock  calculations  is  attributed  to  Von  Neumann 
and  Richtmyer  [6].   They  wanted  to  produce  a  smooth 
transition  across  the  shock  in  all  the  dependent 
variables  when  using  the  partial  differential  equations 
of  motion.   Here,  in  the  same  spirit,  we  introduce  a 
viscous  pressure  ^     wherever  the  pressure  term  p 
appears,  i.e.   substitution  p  +  (j)  ->  p .   In  vector 
notation,  equation  il.2)    then  becomes 


(5.1) 


w,=f   +g  +  r       +s_ 
,t  ,x   ^,y    ,x    ,y 


where  the  new  vectors  are  defined  by 


/?\ 


(5.2) 


r  = 


0 


s  = 


1+' 


The  viscous  pressure  vectors   r  and   s   should  be  very 
small  when  away  from  sharp  compressions,  and  should  attain 
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maximiim  values  in  the  shock  region.  Therefore,  the 
choice  of  the  form  of  (f)  is  dictated  by  the  desire  to 
have  it  small  in  the  smooth  part  of  flow,  and  large  in 
the  transition  or  shock  region.  It  should  also  have  a 
non-zero  contribution  when  the  shocked  flow  is  locally- 
one  dimensional.   A  simple  form  for  the  viscosity  is 


i     _2  /Su  .  dvv 2 


This  form  reduces  to  the  usual  Rlchtmyer  viscosity  in 
one  dimension,  [9]^  and  is  also  invariant  under  a 
rotation  of  the  coordinates.   The  constant  p  has  the 
dimension  of  length. 

We  also  restrict  ^      in  the  following  way: 


^  =   ^'p(f),x-^f^y^^  ^^Cf)^,<0,  (|)^^<0 


(5.3) 


P^p(Cf),,)^ 


P^P((^)  )^ 


=  0 


if  (H)     >  0,  {^)     <  0 

p  ,y  -  '  pSx 


if  (-)   >  0,  (^)   <  0 


if  {^)       >   0,  i^)       >   0 


i.e.,  the  x(y)   component  of  the  viscosity  is  set  equal 
to  zero  when  the  fluid  Is  ixndergoing  a  local  x(y) 
expansion.   As  a  result  of  the  dissipation  terros, 
equation  (5.I)  is  no  longer  a  first  order  quaslllnear 
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partial  differential  equation.  Based  on  previous 
experience  with  the  one  dimensional  form  of  these 
higher  order  equations,  we  expect  the  allowable  time 
Increment  for  stability  for  Its  difference  analog  to  be 
appreciably  smaller  than  required  for  difference  schemes 

(4.1),  (4.2),  and  (4.6). 

2 
The  value  of  p   was  chosen  small  enough  so  that 

the  shock  transition  would  be  as  sharp  as  possible,  but 

large  enough  so  that  oscillations  would  be  damped.   A 

value  of  ,1.9,  seemed  to  offer  the  proper  compromise. 

This  fails  within  the  range  of  values  obtained  by 

Rlchtmyer  [9]  for  the  one  dimensional  case,   p   is 

proportional  to  Ax.   Hence  we  do  not  have  to  consider 

2 

r    and   s    in  the  term  w  ,  ,  — o—  as  we  neglect 
,x        ,-y  ,tt   2 

quantities  of  third  order  and  higher.   The  difference 
scheme  is  then 


(5.M 


w,  .(At)  =  (r   .   (0)-r  .      (0))^ 

+  ^^  .  l^°^-^  .  1  ^°^^^  +  ^i,j^°^ 


subject  to  the  constraints  Imposed  by  (5.3)  and  where 
the  terms  given  on  the  right  hand  side  of  (6.1)  are 
represented  by  f .     .(0). 
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6.  Regular  Net  Points 

A  regular  net  point  of  a  lattice  is  one  which  has 
its  four  nearest  neighbors.   Since  we  are  using  reflec- 
tion principles  on  the  boundary,  a  boundary  point  may 
also  be  a  regular  net  point.   Given  the  initial  data  at 
time  t  =  0,   p(x,y,0),   m(x,y,0),   n(x,y,0),   and 
E(x,y,0),   the  basic  difference  equation  for  the  evalua- 
tion of  the  vector  function  w,   for  the  subsequent  time 
At,   given  equal  net  spacing  Ax  =  Ay  =  A  is 


^,jt^t'  =  «i,jf°'+-i(fi+i,jf°'-fi-i,j(°" 


(6.1)  +-i(s,,j^i(0)-g,_j.,(0)) 

4  2 

+  T~  (-1)^"^  c.(o)  ^     . 

TT  ^  2A'^ 

The  vectors      C.      are   evaluated  respectively  at 

111  1 

(i+-g,j),       (i--2,j),      {1,3+^),      and      (l,j-    ^^i    i.e., 

^1   =   ^i4,j[^fi+l,J-^i,j)    +^fSi,J+l-Si,J-l+gi-M,J+l 

Cp  =   A 

^  1-2. 


=  vi,j(ffi,r^i-i,j)  +^[gi-i,j4-i-gi-i,j-i 

(6.2)  +Si^j+i-gi^j_l]} 

C^      =      B.  .^lii[f.,,  .-f.        ,  .+f  .     ,n  -LT-f-        T  •    LT    ] 

5  l:.J  +  2l^         1+1^  J  1-1.  J  1+1,  J+1         1-1,  J+1' 

-^    fSi,j-hl-Si,j)} 

C,,    =    B.      .     l^-i[f.,,      .    .-f.     .      .    n+f.,n      .-f.     -,      .] 
4  i,j--2[4'-    1+1,  j-1      1-1,  j-1      1+1,  J       1-1,  J 

+  ^ei,rSi,j-i^} 
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The  coefficient  matrices  are  taken  to  be  A.   1  .  = 

1  '- 

-=(A.  ,T  .+A.  .),   etc.   It  is  to  be  noted  that  once 
2^  1+1, J   i,J  ' 

computation  is  completed  at   (i,j),   if  computational 
sequencing  is  by  rows  (columns),   C^  ->  C^  for  point 
(i+l,j),   (Cp  ->  Ch   for  point   (i,j+l)).   This  means 
that,  except  for  regions  of  the  flow  near  the  boundary, 
these  vectors  are  alternately  added  to  and  subtracted 
from  the  flow  and  hence  cancel  In  pairs  in  the  interior 
region.   This  result  follows  as  a  consequence  of 
conservation  differencing,  and  Insures  that  the  mass, 
momentum  and  energy  will  be  conserved  in  the  interior 
region. 

7.   Boundary  Points 

Rigid  walls  are  represented  by  a  string  of  net 
points.   All  such  points  on  these  walls,  which  are 
parallel  to  the  x  or  y  coordinate  axes  are  treated  as 
regular  net  points.   This  is  accomplished  by  using  a 
reflection  principle,  mentioned  earlier,  which  images 
the  flow  properties  onto  a  virtual  line  of  net  points. 
These  points  are  only  used  to  generate  approximations 
to  partial  derivatives  at  the  same  tlm.e  level  for  points 
on  the  neighboring  rigid  wall.   The  conservation 
variables  are  reflected  into  the  line  parallel  to  the 
X  axis  by 

-   20   - 


(except  for  a  necessary  minus  sign  when  reflecting  the 
y  component  of  momentimi) ,  and  parallel  to  the  y  axis  by 

^iTl,j^^^^^^   "  ^i+l,j^^+^^) 

(except  for  a  necessary  minus  sign  when  reflecting  the 
X  component  of  momentum) . 

For  boundaries  which  are  not  parallel  to  either 
coordinates  but  which  make  an  angle  a  with  the  x  axis, 
the  construction  is  somewhat  more  complicated,  although 
the  same  principles  are  used.   In  cases  tested  in  this 
paper,  oblique  surfaces  were  defined  by  a  string  of 
regular  net  points,  i.e.  the  reciprocal  slope  is  an 
integer.   The  slope  of  the  surface  is  known,  and  hence 
allows  the  construction  of  the  local  normal  to  the  surface. 
This  line  segment  is  extended  until  it  cuts  lines  connect- 
ing two  parallel  strings  of  interior  net  points  at   L 
and  L    (see  figure) .   Quadratic  interpolation  of  the 
conservation  variables  at  these  two  points  define,  with 
with  the  help  of  the  original  (or  quadratically  interpolated) 
boundary  point,  on  the  surface  L  ,   a  quadratic  function 
along  the  normal  to  the  surface.   In  order  to  compute  this 
quadratic  function,  one  must  evaluate  a  first  derivative 
to  second  order  accuracy,  and  a  second  derivative,  having 
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to  use  unequal  spacing  along  the  normal  L   to  the 
surface.   Then  these  derivatives  are  given  in  terms  of 
the  eccentricity   5~  =  (L  -  L  )/(L  -  L__)  . 


^2 

0     W 


L         +         o      - 
o 


^   (w^  -  (l+5)w^  -  6Wt-  )     ^ 


L         ^+        ^o      -  dib+DAL"^ 
o 


It  was  found  that  these  formulas  did  quite  well  for 
most  image  points  in  the  wedge.   Several  points,  however, 
gave  rise  to  values  of   5~  »  1,   i.e.,   L  %  L_ .   After 
several  tens  of  cycles  of  computation,  instabilities 
(negative  densities)  were  noted.   At  these  points  then, 
rather  than  interpolate  the  functional  values  at  the 
first  two  lines  of  net  points,  the  second  and  third  line 
of  net  points  were  used.   This  appeared  to  solve  the 
problem,  since  the  value  of  the  eccentricity  was   6  %  1. 
This  phenomenon  is  apparently  similar  to  that  observed  by 
John  Gary  in  his  shock  fitting  work  [10].   As 
Robert  Rlchtmyer  pointed  out  [7],  whenever  two  points  in 
a  computation  come  much  closer  than  a  characteristic 
length  equal  to   cAt,   instability  may  be  the  result. _ 
For  the  momentum  variables 


22   - 


in(t+At)     =  j'2(ni(t+At)cosa+ n(t+At)sina)cosa-  m(t+At)j"  L* 

image  point  ^  ' 

n(t4-At)     =  |2(m(t+At)cosa+ n(t+At)sina)sina- n(t+At)^L* 

Image  point  ^  ' 


The  density  p   and  energy  E  may  be  reflected  by  the  rules: 


p(t+At)   =  p(t+At) 
image  point 


E(t+At)   =   E(t+At) 


Finally,  we  consider  the  downstream  boundary.   Since 

an  x-difference  at  the  downstream  boundary   (i=i   ) 

"^      max 

cannot  be  evaluated  in  the  usual  manner,  we  use  extrapola- 
tion formulas.  Using  data  at  points   i   ,  i    -,  ,  and 

max   max-1 

ir.,0^  o  (along  the  rows   j-l,J,j+l),  we  may  form 
md.x—  ^ 

x-differences  of  the  vectors   f  and  g.   The  extrapolated 

■X-     * 

values  are  denoted  by  f  ,  g  ,  etc.   Then  at   i=i  ^  ,  the 

■^  x'    ^x'  max' 

downstream  boundary  is  computed  by 


4 


w(At)  =  w(o)  +  At(f%g^)  +^y~{-i)^-'^c](^) 


y 


^=1 


The  vectors  C„      contain  extrapolated  values  of  f  and  g 

i/  XX 

wherever  these  differences  are  needed  in  (6.2). 
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8.   Shock  Relations 

The  Integral  form  of  the  conservation  equations  are 
to  be  satisfied  across  the  shock.   These  equations,  after 
some  algebraic  manipulations,  may  be  put  into  polynomial 
form,  i.e.,  a  cubic  equation  for  the  local  tangent  of  the 
shock  angle 

tan-^C  -  C„tan  c  +  C-,  tancf  +  C   =0 
2  1         o 

where 

C^     =      2(M^-l)/tana(2+(7-l)M^) 

0^     =    (2+(7+l)M^)/(2+(7-l)M^) 

Cq  =  2/tana(2+(7-l)M^) 

The  angle  the  flow  makes  with  the  horizontal  is  a, 

while  the  Mach  number  approaching  the  shock  is  M  >  1. 

It  follows  from  Descartes  rule  of  signs  that  this  equation 

has  no  more  than  two  positive  roots,  and  one  negative  root. 

Since  the  discriminant  is  negative  for  certain  values  of 

a   and  M,  we  expect  a  strong   (7  ,   weak  o"  ,   and  an 
'  ^  "^    s  w 

impossible  solution  c.^,  {a     >   <5  ,      cTt  ^  ^o)«   The 

_L        S      W      J-      S 

solution  cC   will  produce  a  density  and  pressure  decrease 
across  the  shock,  a  result  that  leads  to  a  local  reduction 
of  entropy  with  no  corresponding  increase  elsewhere  in  the 
universe.   Hence,  a  contradiction  of  the  second  law.   The 
choice  of  the  remaining  two  solutions  is  dictated  by  the 
boundary  condition  on  the  pressure  downstream  of  the  shock, 
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Then  if  the  subscript  1  (2)  denotes  the  position  before 
(after)  the  discontinuity,  the  density,  velocity  and 
energy  ratios  are  solved  for,  viz. 


Pi 

(tanC-tana) 

1 

P2 

1  +  tana'  tana 

tano- 

^1 

cos(o-  -a) 

^2 
^1 

coscr 

P2,7+1        P2^ 
Pl^7-1        Pi^ 

^2 

7+1      P  2 
7-1      Pi 

These  equations  were  used  to  generate  the  exact  solution 
of  the  flow  field  for  the  regular  oblique  reflection 
patterns  and  these  were  compared  to  the  numerical 
solution  (see  Table  2) .   The  shock  angles  defined  the 
extent  of  each  plecewlse  constant  region  while  the  jump 
conditions  gave  the  value  of  the  corresponding  state 
properties.   However,  when  M  is  too  small  for  a  given 
a,  or  a   too  large  for  a  given  M,   so  as  to  preclude  a 
regular  reflection  solution,  a  Mach  reflection  occurs. 
In  this  case,  the  above  equations  give  complex  roots  and 
cannot  be  used  to  predict  the  flow  pattern  for  the 
resulting  triple  shock  wave  (sometimes  called  lambda 
shock)  configuration.   As  a  result,  the  Mach  reflection 
which  was  also  generated  as  a  time  dependent  problem, 
had  no  exact  solution  for  comparison.   Given  some 
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arbitrary  initial  data,  but  using  correct  boundary 
conditions  at  x  =  -.1  (-00),  i.e. 

M_^   <  1.87  for  7  =  1.4  and  a  =  tan~^(l/5); 

the  result  would  be  a  numerical  solution  with  three 
intersecting  shocks:   Incident  or  primary  shock, 
reflected  shock,  and  the  Mach  stem.   These  asymptotic 
solutions  are  shown  in  the  figures  for  various  values 
of  M  (=  1.85,  1.7,  1.6  and  I.5  with  7  =  1.4).   The 
complete  time  dependent  solution  for  M   =  1.6  and 
7  =  1,2   is  also  shown  to  indicate  the  transient  nature 
of  the  Mach  reflection  calculation. 

Although  the  above  relations  are  convenient  for  the 
present  problem,  the  general  form  of  the  Jump  conditions 
may  be  obtained  from  the  differential  equations  (2.2). 
If  the  left-hand  side  is  integrated  over  a  volume  region 
of  space,  while  the  right-hand  side  is  integrated  over 
the  surface  defining  this  volume  region,  the  result 

S[w^]   =   [F^  •  n] 

is  the  Rankine-Hugoniot  relations.   Here,   n  Is  the  unit 
normal  vector  to  the  discontinuity,   P.   is  the  vector 
(f-^g-)j   and  S   is  the  normal  speed  of  the  discontinuity. 
[a]  indicates  the  jump  in  A.   Every  weak  solution  must 
satisfy  these  j\imp  conditions  across  a  line  of  discontinuity, 
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This  fact  was  used  to  test  the  two  dimensional  niomerlcal 
method  as  far  as  Its  ability  to  reproduce  one  dimensional 
shocks.   If  the  correct  initial  conditions  were  given, 
i.e.  a  stationary  one  dimensional  shock  solution,  the 
difference  equations  were  identically  satisfied  for  all 
t   if  at  t  =  0,   the  flow  field  has  values 

-00  <  X  <  0        w  =  w(x<0,0)   =  w_(0) 

0  _<  X  <  +00      w  =  w(x>O,O0,  =   w  (0) 

To  see  this,  the  Rankine-Hugoniot  condition  in  this  case 
(S=0)   is 

(8.1)  [f]   =  0 

i.e.,   f  _^  =  f_. 

The  two  dimensional  equations  reduce  to  the  one 
dimensional  form 

2 

w(x,t+At)   =   w(x,t)  +  At  f/s  + -^  (Af_) 

X    d  X  X 

In  front  of  and  behind  the  shock,  the  state  is  constant. 
Therefore,  all  derivatives  in  the  second  order  term 
vanish.   Across  the  shock  (8.1)  is  satisfied  so  that  at 
the  shock  a  'centered  difference  in   f,   which  is 
numerically  equal  to  the  jiimp  in   f,   results  in  the 
vanishing  of  the  difference  appearing  in  the  first  order 
term. 
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9.   Initial  Conditions 

The  four  difference  schemes  were  tested  on  a  problem 
of  ordinary  shock  reflection.   The  flow  is  supersonic  and 
plane  until  a  constriction  or  wedge  appears  in  the 
channel.   In  this  region,  an  incident  shock  changes  the 
flow  direction  through  an  angle   a,   so  that  it  is 
parallel  to  the  wedge.   A  secondary,  or  reflecting  shock, 
changes  the  flow  direction  by  an  equal  amount  so  that  it 
is  again  parallel  to  the  upper  wall.   Since  these  initial 
conditions  represent  an  exact  solution,  it  was  felt  that 
a  comparison  of  the  difference  solution  of  this  steady 
flow  problem  with  the  exact  solution  would  be  a  good 
measure  of  the  computing  methods . 

An  additional  problem  dealing  with  the  time 
dependent  approach  of  the  flow  to  steady  state  ordinary 
reflection  solution  was  performed  with  arbitrary  initial 
data  so  as  to  ascertain  the  rate  of  convergence  of  the 
method.   The  incident  shock  was  severely  perturbed  as  it 
was  set  standing  vertically  (cr  =   ^)   at  the  toe  of  the 
wedge   (a  was  given  the  value  tan"  (I/5)).   In  addition, 
the  jump  conditions  across  this  shock  were  made  to  differ 
from  those  predicted  by  the  Hugoniot  conditions.   The 
reflected  shock  was  completely  eliminated  In  the 
initialization.   Here,  the  initial  conditions  were 
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-00  <  X  j<  0 
and 
0  <  y  <  y  ,  T 

00  >  X  >  0 
and 

"'wedge-'^  —  •^wall 


w  =   w(x  <  0,y,0) 


w  (0) 


J 


w  =   w(x  >  0,y,0)   =   w  (0) 


The  results  are  shown  in  the  accompanying  figures,  for 
values  of  the  Mach  number  M   =2  and  2.2,  where  the  density 

-00  '  "^ 

and  pressure  fields  are  plotted.   An  indication  of  the  rate 

of  convergence  of  the  time  dependent  solution  is  shown  for 

M  _=  2.2  case. 
-oo 


10.   Numerical  Tests  of  Stability 

For  all  the  methods  tested,  the  C.F.L.  stability 
number   (At/Ax)^^^^^^  =  ^^^/^y^C.F.L.  "  l/(|u|+c)  =  0.66?. 
The  difference  scheme  associated  with  the  triple 
viscosity  (4.6)  should  exhibit  stability  at  this  value 
of  the  stability  number.   A  series  of  tests  were  run  with 
this  difference  scheme  being  applied  to  the  ordinary 
reflection  problem  with  varying  stability  numbers.   It 
■was  found  that  at  .97 (At/Ax) „  .^  ^      instability  occurred 
at  29  cycles;  at  .86  (At|Ax)p  „  ^  instability  occurred 
at  87  cycles,  but  at  .80  (At/Ax)„  „  ^   the  problem  was 
run  to  400  cycles  with  no  osciallations  appearing  in  the 
flow.   The  term  instability,  as  it  is  used  here,  means  the 
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appearance,  at  one  or  more  net  points,  of  negative  values 
of  the  density.   The  probable  reason  for  the  difference  of 
the  theoretical  value  compared  with  the  empirical  value 
lies  In  two  basic  assumptions  of  the  stability  analysis; 
a)  The  flow  Is  smooth  so  that  the  coefficient  matrices 
A  and  B  are  practically  constant;  b)  The  boundaries 
have  no  effect  on  the  form  of  the  amplification  matrix. 
As  far  as  a)  Is  concerned  the  m.atrlces  A  and  B  are  not 
even  remotely  constant  in  the  shock  transitions.   The 
higher  derivatives,  associated  with  the  triple  viscosity, 
take  on  large  values  in  this  region.   A  test  was  made  in 
which  the  two  third  order  terms  were  set  equal  to  zero, 
while  retaining  only  the  fourth  order  term.   At 
.86  (At/Ax)„  -I-,  -r  ,    the  computation  appeared  to  be 
stable  to  120  cycles.   No  further  tests  were  run  with 
this  truncated  scheme.   As  far  as  b)  is  concerned,  little 
is  known  about  the  behavior  of  the  boundaries  on  the 
stability  properties  of  these  difference  schemes. 

Equation  (3.3)  was  also  run  at  varying  fractional 
values  of  the  Courant  ni;iinber.   In  this  scheme,  the  highest 
order  difference  is  two,  as  contrasted  with  four  of  the 
previous  case.   At  .86  (At/Ax)^   ^  =  .86  (At/Ay)^  ^  ^  > 
the  computation  was  carried  out  to  400  cycles  with  no 
instability.   This  was  a  time  interval  greater  than  that 
prediced  in  [4]  by  a  factor  of  approximately  2. 
Instability  resulted  if  larger  time  intervals  were  used. 
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It  is  clear  that  in  the  stability  proof  of  [k]    the  use 
of  inequality  estimates  leads  to  conservative  results. 
Equations  (30)  and  iJ>A)    were  also  found  stable  at 
.86(At/Z\x)„  -r,  T  >  ^'^t  unstable  in  20  cycles  at 
.9^  (At/Ax)„  „  T-    The  last  difference  method  tested 
(5.4),  however,  was  found  unstable  at  all  these  large 
values  of  the  time  step.   Instead,  to  insure  the  stability 
of  the  scheme,  increments  less  than  0.1  (At/Ax)„  „  ^   were 
necessary.   A  comparison  of  these  methods  is  shown  in  a 
test  case  in  Table  1. 

11 .   NiAmerical  Results 

The  numerical  scheme  (6.1)  was  tested  for  varying 
values  of  the  Mach  number  of  the  fluid  entering  the 
channel,  viz.   1-3  jf  '^-nn—  5*^'   ^°^  ^^®  wedge  angle 
used,  the  flow  followed  one  of  two  solutions  which  are 
distinct  and  which  depend  quite  strongly  on  the  value 

of  the  Mach  n\imber,  i.e.  for  M   <  I.87,  Mach  reflec- 

'  -00—     ' 

tions  occurred,  while  for  larger  values  of  the  Mach 
niimber,  the  asymptotic  solution  was  that  of  ordinary 
reflection.   The  methods  used  were  able  to  predict 
quite  sharply  this  limiting  Mach  number.   Table  2  gives 
the  results  of  the  shock  angles  for  both  sets  of 
solutions.   Comparison  of  the  computed  data  with  the 
exact  solution  for  ordinary  reflection  seems  to  be  good. 
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The  Mach  reflection  calculation  was  not  compared  to 
experimental  data  but  the  behavior  seems  to  be  quali- 
tatively correct.   A  steady  state  solution  for  values 
of  M   <  1.4   could  not  be  obtained.   The  reason  seems 

-00  — 

to  be  associated  with  a  combination  of  factors.   The 
rapidly  shrinking  triangular  region  under  the  triple 
point  combined  with  the  technique  of  computing  deriva- 
tives at  the  surface  of  the  wedge  leads  to  inaccuracies. 
The  reflection  method  that  is  used  on  the  wedge  is  a 
second  order  scheme  and  hence  will  be  inaccurate  when 
there  are  large  variations  in  the  flow  properties,  i.e., 
when  many  shocks  are  present  in  a  small  region.   In 
addition,  for  these  low  values  of  Mach  number,  the  long 
Mach  shock  means  that  the  flow  is  basically  subsonic  in 
the  downstream  region,  a  fact  which  allows  disturbances 
(which  might  be  generated  numerically  by  the  extrapola- 
ting procedure)  to  propagate  upstream  from  that  boundary. 

In  addition,  the  sllpline  emerging  from  the  triple 
point,  and  which  divides  the  flow  downstream  in  the  Mach 
reflection,  could  not  be  numerically  resolved  in  these 
calculations.   This  was  probably  due  to  the  weakness  of 
the  shocks  in  this  study,  i.e.,  the  strength  of  the 
sllpline  depends  upon  the  difference  of  entropy  across 
this  line,  and  the  entropy  production  through  a  weak 
shock  behaves  as  the  cube  of  the  shock  strength. 
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TABLE  1 
Comparison  of  the  Four  Difference  Methods  of  the  Density 


at  Upper 

Boundary  of 

Oblique  Reflection  M_^  = 

:  yu^+v^/c  = 

Equation 

Equation 
(3.3)+(3.4) 

Equation 
(4.6) 

Equation 
(5.4) 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

exact 

density 

1.0 

1.000 
1.000 
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1.000 
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1.000 
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.999 
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.998 

.998 

.999 

1.002 
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.997 

1.000 

1.007 

1.005 

1.011 

1.000 

.996 

.993 

.988 

.987 

.982 

.970 

.986 

.991 

.964 

.948 

.957 

1.068 

.962 

.975 

.995 

1.239 

exact 

1.246 

1.283 

1.199 

1.499 

shock  — 
position 

""^.072 
2.509 

2.035 

2.479 

2.097 
2.678 

1.796 

2.042 

exact 

2.319 

2.356 

2.248 

2.214 

density 

2.296 

2.308 

2.265 

2.323 

2.277 

2.327 

2.320 

2.338 

2.342 

2.318 

2.314 

2.319 

2.306 

2.319 

2.317 

2.306 

2.246 

2.337 

2.331 

2.347 

2.216 

2.334 

2.330 

2.366 

2.220 

2.341 

2.337 

2.369 

2.212 

TABLE  2 
Results  of  Oblique  Reflection 

(7  =  1.^) 


M  Wedge  Angle 

—00     °    ° 


2.0 
2.2 
5.0 


tan"^  (h 


Incident  Shock 
Computed 

40.2° 

^0 


Reflected  Shock 


37.2^ 
27.5 


o 


Exact 
40.8° 
37.2° 
28.6° 


Computed 

45.8° 

47.0° 

not 
computed 


Exact 
54.0° 
46.6° 
55.8° 


M     Wedge  Angle 
-00     °    ° 

1.85    tan"-^  (-i) 

1.7 
1.6 

1.5        " 
1.6(7=1.2)   " 


Mach  Reflection 
(7  -   1.4) 

Mach  Stem  0^ 


90.0 

90.0 

90.0^ 

90.0^ 

90.0 


o 


Incident 
Shock  0^ 


o 


46.2 

48.7^ 

54.6' 

59.0' 

52.5^ 


o 


Reflected 
Shock  9^ 

ee""  -  80° 

61°- 73° 

65.3° 


70.2 
72.4^ 


o 


II  M  M  II  I' 


List  of  Figures 

l.a,b   Asymptotic  solution  for  ordinary  reflection  for 

M_^=  2.0,  Density  field  and  pressure  field. 

Overall  density  ratio  =  2.3,  exact  =  2.27 

(At/A  =  Ahj>) 
l.c     Approach  to  the  asymptotic  solution  for  ordinary 

reflection  for  M_  =  2.0  for  the  momenta  m  and  n 

near  the  upper  wall  (y  =  O.75)   At/A  =  0.443. 
2. a- J   Approach  to  the  asymptotic  solution  for  ordinary 

reflection  for  M_^=  2.2.   Density  field  and 

pressure  field  (Pig.  2j).   Overall  average 

density  ratio  =  2.33,  exact  =  2.34.   Figures 

correspond  to  the  density  field  at  35,  70,  IO5, 

140,  210,  350,  420,  560,  700  times  At. 

(At/A  =  .233) 
3.a,b   Asymptotic  solution  for  Mach  reflection  for 

M_  =  1.85.   Density  field  and  pressure  field. 

Figures  correspond  to  5OO  At   (At/A  =  .48) 
4.a,b   Asymptotic  solution  for  Mach  reflection  for 

M_  =  1.7.   Density  field  and  pressure  field. 

Figures  correspond  to  420  At   (At/A  =  .229). 
5.a,b   Asymptotic  solution  for  Mach  reflection  for 

M_  =  1.6.   Density  field  and  pressure  field. 

Figures  correspond  to  420  At   (At/A  =  .233). 
6.a-h   Approach  to  the  asymptotic  solution  for  Mach 

reflection  for  M   =  1.6  with  7  =  1.2. 

Density  field  and  pressure  field  (Fig.  6h) . 

Figures  correspond  to  the  density  field  at 

70,  140,  210,  350,  420,  560,  700  times  At 

(At/A  =  .229) . 
7.a,b   Asymptotic  solution  for  Mach  reflection  for 

M_  =  1.5.   Density  field  and  pressure  field. 

Figures  correspond  to  490  At   (At/A  =  .233). 
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